#ifndef RIGID_HPP
#define RIGID_HPP
#include "cell.h"
#include "utils.hpp"
template <int NATOMS, int... IJS>
struct rigid_type {
  static constexpr int ijs[sizeof...(IJS)] = {IJS...};
  static constexpr int ncons = sizeof...(IJS) / 2;
  static constexpr int natoms = NATOMS;
  static constexpr auto iaseq = utils::make_iseq<int, natoms>();
  static constexpr auto icseq = utils::make_iseq<int, ncons>();
  template <int side, int... ICs>
  static constexpr auto make_ij(utils::iseq<int, ICs...> icseq) {
    return utils::iseq<int, ijs[ICs * 2 + side]...>{};
  }
  static constexpr auto iseq = make_ij<0>(icseq);
  static constexpr auto jseq = make_ij<1>(icseq);
  static constexpr auto ic2seq = utils::make_iseq<int, ncons * ncons>();
  static constexpr auto ic3seq = utils::make_iseq<int, ncons * ncons * ncons>();
  template <int... ICsq>
  static constexpr auto make_icmat(utils::iseq<int, ICsq...> icsqseq) {
    return utils::iseq<int, ICsq / ncons...>{};
  }
  template <int... ICsq>
  static constexpr auto make_jcmat(utils::iseq<int, ICsq...> icsqseq) {
    return utils::iseq<int, ICsq % ncons...>{};
  }
  static constexpr auto icmatseq = make_icmat(ic2seq);
  static constexpr auto jcmatseq = make_jcmat(ic2seq);
  static constexpr real sign(int ic, int jc) {
    int i0 = ijs[ic * 2];
    int j0 = ijs[ic * 2 + 1];
    int i1 = ijs[jc * 2];
    int j1 = ijs[jc * 2 + 1];
    if (i0 == i1 || j0 == j1)
      return -1;
    else
      return 1;
  }

  static constexpr int conn(int ic, int jc) {
    int i0 = ijs[ic * 2];
    int j0 = ijs[ic * 2 + 1];
    int i1 = ijs[jc * 2];
    int j1 = ijs[jc * 2 + 1];
    if (i0 == i1)
      return i0;
    if (i0 == j1)
      return i0;
    if (j0 == i1)
      return j0;
    if (j0 == j1)
      return j0;
  }
};

using rigid2 = rigid_type<2, 0, 1>;
using rigid3 = rigid_type<3, 0, 1, 0, 2>;
using rigid4 = rigid_type<4, 0, 1, 0, 2, 0, 3>;
using rigid3angle = rigid_type<3, 0, 1, 0, 2, 1, 2>;
enum rigid_code {
  RIGID2 = 1,
  RIGID3 = 2,
  RIGID4 = 3,
  RIGID3ANGLE = 4
};
#include "memory_cpp.hpp"
#include "elemtab.h"
struct rigid_type_table {
  int *base, *types, *bonds, *angles, ntypes;

  rigid_type_table(int ntypes) {
    this->ntypes = ntypes;
    this->base = esmd::allocate<int>(ntypes + ntypes * ntypes + ntypes * ntypes * ntypes, "rigid/type table");
    this->types = this->base;
    this->bonds = this->types + ntypes;
    this->angles = this->types + ntypes * ntypes;
    for (int i = 0; i < ((ntypes + 1) * ntypes + 1) * ntypes; i ++){
      this->base[i] = 0;
    }
  }
  rigid_type_table(const char *condstr, elemtab_t *elemtab) : rigid_type_table(elemtab->cnt) {
    parse(condstr, elemtab);
  }
  bool is_bond_rigid(int ti, int tj) {
    if (this->bonds[ti * ntypes + tj] == -1 || this->bonds[tj * ntypes + ti] == -1) return false;
    return this->types[ti] || this->types[tj] || this->bonds[ti * ntypes + tj] || this->bonds[tj * ntypes + ti];
  }
  bool is_angle_rigid(int ti, int tj, int tk) {
    return this->angles[(ti * ntypes + tj) * ntypes + tk];
  }
  void parse(const char *condstr, elemtab_t *elemtab);
  void find_clusters(cellgrid_t *grid, real *bond_len, real *angle_len);
};
#include "comm.h"
template <class T>
void rigid_coordinates(cellgrid_t *grid, mpp_t *mpp, T &rigid);
template <class T> void rigid_force (cellgrid_t *grid, real dtfsq, T &rigid);
template <class T>
void rigid_force_listed(cellgrid_t *grid, real dtfsq, T &rigid);
template <typename T>
void rigid_setup(cellgrid_t *grid, mpp_t *mpp, real dt, real ftm2v, T &rigid);
template <typename T> void rigid_post_force (cellgrid_t *grid, mpp_t *mpp, real dt, real ftm2v, T &rigid);
template <typename T> void rigid_post_force_listed (cellgrid_t *grid, mpp_t *mpp, real dt, real ftm2v, T &rigid);

#endif
